#like 2.1, but with separate tables for lead and no-lead
#   donor cases, ad geographic specification

#separate lead donor and no-lead donor findings
#ad geographic specification

#produce figures
rm(list=ls())
library(epicalc)
library(MASS)
zap()

source("C:/Users/Martin Steinwand/Documents/R/MyRoutines/LatexTable.R")
setwd("C:/Users/Martin Steinwand/Documents/Project Lead Donor/Analysis/R&R.IO")
nm <- 15 #number of models

#create list to hold results
res <- vector("list", nm)
names(res) <- c("AllNoLDWev", "AllStrongLDWev", "AllStrongLDWld","NGONoLDWev",
                "NGOStrongLDWev", "NGOStrongLDWld", "GovNoLDWev", "GovStrongLDWev", "
                GovStrongLDWld")
fnames <- c("AllPop.pooled.wev","AllPop.strongld.2.2.wev", 
            "AllPop.strongld.2.2.wld",
            "NGOPop.pooled.wev","NGOPop.strongld.2.2.wev", "NGOPop.strongld.2.2.wld",
            "GovPop.pooled.wev","GovAidpop.strongld.2.2.wev", "GovAidpop.strongld.2.2.wld",
            "AllPop.pooled.geo", "AllPop.strongld.2.2.geo", "NGOPop.pooled.geo", 
            "NGOPop.strongld.2.2.geo", "GovPop.pooled.geo", "GovAidpop.strongld.2.2.geo")
for(i in 1:length(fnames)){
    load(fnames[i])
    res[[i]] <- get(objects()[5])
    rm(list=objects()[5])
}

#obtain draws from sampling distribution for all models

reps <- 10000
#reps <- 1000

getbds <- function(no){
    pars.mc <- matrix(NA, reps, length(res[[no]]$par))    
    for (i in 1:reps){
            pars.mc[i,] <- mvrnorm(1, res[[no]]$par, res[[no]]$vacov)
    }
    pars.out <- apply(pars.mc, 2, quantile, probs=c(.025, .5, .975))
    pars.out[2,] <- res[[no]]$par 
    return(pars.out)
}

bds <- lapply(1:nm, getbds)


##############################
## no lead donor

modname <- c("All channels",
                "NGO aid",
                "Government aid",

               "All channels",
                "NGO aid",
                "Government aid")


    #note: models 1, 4, 7, 10, 12, 14 have 17 parameters, the others 18. Solution: suppress
    # parameter on position 4 (donor dummies) in graphical representation
modord <- c(1, 4, 7, 10, 12, 14)
plotout <- array(NA, c(3, ncol(bds[[1]]), length(modname)))
for (j in 1:length(modname)){
    if(ncol(bds[[modord[j]]])>ncol(bds[[1]])){
        plotout[,1:3,j] <- unlist(bds[modord[j]])[1:9]
        plotout[,4:ncol(bds[[1]]),j] <- unlist(bds[modord[j]])[13:length(unlist(bds[modord[j]]))]
    } else plotout[,,j] <- unlist(bds[modord[j]])
}

#no ld: geo vs even 
#       All -- even wins
#       NGO -- geo wins
#       Gov -- even wins
#ld: ld vs geo 
#       All -- ld wins
#       NGO -- geo wins 
#       Gov -- ld wins
#    ld vs even -- as before
#       All -- ld wins
#       NGO -- even wins
#       Gov -- ld wins
#    even vs geo 
#       All -- even wins
#       NGO -- even wins
#       Gov -- even wins
#

#pick clarke test winners
winord <- c(1, NA, 3, NA, 5, NA) 
winout <- array(NA, dim(plotout))
for (j in 1:(length(modname))){
    winout[,,j] <- plotout[,,winord[j]]
}

###
#rho 
y.axis <- c(7:5, 3:1)

low <- min(plotout[,1,], na.rm = TRUE)
hgh <- max(plotout[,1,], na.rm = TRUE)

par(mar=c(2, 11, 2, .5), xpd=TRUE) 
plot(c(low, hgh), c(1,max(y.axis)), type = "n", 
    axes=F,
    xlab = "", 
    ylab = "", #turn off axes and labels. 
        xlim=c(low, hgh), c(1, max(y.axis)), xaxs = "r",
         main = expression(rho)) 
#the 3 lines below create horiztonal lines for 95% confidence intervals, and vertical ticks for 90% intervals
segments(plotout[1,1,], y.axis, plotout[3,1,], y.axis, 
            lwd =  1.5, lty=2)#coef +/-1.96*se = 95% interval, lwd adjusts line thickness
points(plotout[2,1,], y.axis, pch = 19, cex=.6)
#emphasize winners
segments(winout[1,1,], y.axis, winout[3,1,], y.axis, lwd =  2)#coef +/-1.96*se = 95% interval, lwd adjusts line thickness
points(winout[2,1,], y.axis, pch = 19, cex=.9)
#axis for drawing axis
axis(1, at = c(low,hgh), tick=T, labels=c("",""),lwd.ticks=0, cex.axis = .8, mgp = c(2,.5,0))
axis(1, at = c(-.2, 0, .2, .4), tick = T, cex.axis = .8, mgp = c(2,.5,0))
axis(2, at = y.axis, label = modname, las = 1, tick = T, cex.axis = .8) 
lines (c(0,0), c(0,max(y.axis)+1), lwd=1, lty=3)
legend(-.86, 11.2,c("Even weights"), bty="n")
legend(-.86, 7.2,c("Even Weights"), bty="n")
legend(-.86, 3.2, c("Geographic weights"), bty="n")
legend(.25, 3, c("Best fitting models", "Loser in Clarke test"), lty=c(1,2), lwd=1)
#Still to do: highlight which of the models wins 


#save
postscript(file="RhoNoLD.R&R.eps",onefile=FALSE,
 horizontal=FALSE, paper="letter", height=6.1, width=8)
par(mar=c(2, 10, 2, 2), xpd=TRUE) 
plot(c(low, hgh), c(1,max(y.axis)), type = "n", 
    axes=F,
    xlab = "", 
    ylab = "", #turn off axes and labels. 
        xlim=c(low, hgh), c(1, max(y.axis)), xaxs = "r",
         main = "") 
#the 3 lines below create horiztonal lines for 95% confidence intervals, and vertical ticks for 90% intervals
segments(plotout[1,1,], y.axis, plotout[3,1,], y.axis, 
            lwd =  1.5, lty=2)#coef +/-1.96*se = 95% interval, lwd adjusts line thickness
points(plotout[2,1,], y.axis, pch = 19, cex=.6)
#emphsize winners
segments(winout[1,1,], y.axis, winout[3,1,], y.axis, lwd =  1.5)#coef +/-1.96*se = 95% interval, lwd adjusts line thickness
points(winout[2,1,], y.axis, pch = 19, cex=.7)
#axis for drawing axis
axis(1, at = c(low,hgh), tick=T, labels=c("",""),lwd.ticks=0, cex.axis = .8, mgp = c(2,.5,0))
axis(1, at = c( -.1, 0, .1, .2, .3), tick = T, cex.axis = .8, mgp = c(2,.5,0))
axis(2, at = y.axis, label = modname, las = 1, tick = T, cex.axis = .9) 
lines (c(0,0), c(0,max(y.axis)+.2), lwd=1, lty=3)
legend(-.25, 7.7,c("Even Weights"), bty="n", text.font=2)
legend(-.25, 3.7, c("Geographic weights"), bty="n", text.font=2)
legend(.15, 1.78, c("Best fitting models", "Loser in Clarke test"), lty=c(1,2), lwd=1,
            cex=.8)
dev.off()




##############################
## with lead donor

modname <- c("All channels",
                "NGO aid",
                "Government aid",

               "All channels",
                "NGO aid",
                "Government aid",

               "All channels",
                "NGO aid",
                "Government aid")

fnames <- c("AllPop.pooled.wev","AllPop.strongld.2.2.wev", 
            "AllPop.strongld.2.2.wld",
            "NGOPop.pooled.wev","NGOPop.strongld.2.2.wev", "NGOPop.strongld.2.2.wld",
            "GovPop.pooled.wev","GovAidpop.strongld.2.2.wev", "GovAidpop.strongld.2.2.wld",
            "AllPop.pooled.geo", "AllPop.strongld.2.2.geo", "NGOPop.pooled.geo", 
            "NGOPop.strongld.2.2.geo", "GovPop.pooled.geo", "GovAidpop.strongld.2.2.geo")

    #note: models 1, 4, 7, 10, 12, 14 have 17 parameters, the others 18. Solution: suppress
    # parameter on position 4 (donor dummies) in graphical representation
modord <- c(3, 6, 9, 2, 5, 8, 11, 13, 15)
plotout <- array(NA, c(3, ncol(bds[[1]]), length(modname)))
for (j in 1:length(modname)){
    if(ncol(bds[[modord[j]]])>ncol(bds[[1]])){
        plotout[,1:3,j] <- unlist(bds[modord[j]])[1:9]
        plotout[,4:ncol(bds[[1]]),j] <- unlist(bds[modord[j]])[13:length(unlist(bds[modord[j]]))]
    } else plotout[,,j] <- unlist(bds[modord[j]])
}

#no ld: geo vs even 
#       All -- even wins
#       NGO -- geo wins
#       Gov -- even wins
#ld: ld vs geo 
#       All -- ld wins
#       NGO -- geo wins 
#       Gov -- ld wins
#    ld vs even -- as before
#       All -- ld wins
#       NGO -- even wins
#       Gov -- ld wins
#    even vs geo 
#       All -- even wins
#       NGO -- even wins
#       Gov -- even wins
#

#pick clarke test winners
winord <- c(1, NA, 3, NA, 5, NA, NA, NA, NA) 
winout <- array(NA, dim(plotout))
for (j in 1:(length(modname))){
    winout[,,j] <- plotout[,,winord[j]]
}

###
#rho 
y.axis <- c(11:9, 7:5, 3:1)

low <- min(plotout[,1,], na.rm = TRUE)
hgh <- max(plotout[,1,], na.rm = TRUE)

par(mar=c(2, 11, 2, .5), xpd=TRUE) 
plot(c(low, hgh), c(1,max(y.axis)), type = "n", 
    axes=F,
    xlab = "", 
    ylab = "", #turn off axes and labels. 
        xlim=c(low, hgh), c(1, max(y.axis)), xaxs = "r",
         main = expression(rho)) 
#the 3 lines below create horiztonal lines for 95% confidence intervals, and vertical ticks for 90% intervals
segments(plotout[1,1,], y.axis, plotout[3,1,], y.axis, 
            lwd =  1.5, lty=2)#coef +/-1.96*se = 95% interval, lwd adjusts line thickness
points(plotout[2,1,], y.axis, pch = 19, cex=.6)
#emphasize winners
segments(winout[1,1,], y.axis, winout[3,1,], y.axis, lwd =  2)#coef +/-1.96*se = 95% interval, lwd adjusts line thickness
points(winout[2,1,], y.axis, pch = 19, cex=.9)
#axis for drawing axis
axis(1, at = c(low,hgh), tick=T, labels=c("",""),lwd.ticks=0, cex.axis = .8, mgp = c(2,.5,0))
axis(1, at = c(-.2, 0, .2, .4), tick = T, cex.axis = .8, mgp = c(2,.5,0))
axis(2, at = y.axis, label = modname, las = 1, tick = T, cex.axis = .8) 
lines (c(0,0), c(0,max(y.axis)+1), lwd=1, lty=3)
legend(-.86, 11.2,c("Lead donor weights"), bty="n")
legend(-.86, 7.2,c("Even weights"), bty="n")
legend(-.86, 3.2, c("Geographic weights"), bty="n")
legend(.25, 3, c("Best fitting models", "Loser in Clarke test"), lty=c(1,2), lwd=1)
#Still to do: highlight which of the models wins 


#save
postscript(file="RhoLD.R&R.eps",onefile=FALSE,
 horizontal=FALSE, paper="letter", height=6.1, width=8)
par(mar=c(2, 10, 2, 2), xpd=TRUE) 
plot(c(low, hgh), c(1,max(y.axis)), type = "n", 
    axes=F,
    xlab = "", 
    ylab = "", #turn off axes and labels. 
        xlim=c(low, hgh), c(1, max(y.axis)), xaxs = "r",
         main = "") 
#the 3 lines below create horiztonal lines for 95% confidence intervals, and vertical ticks for 90% intervals
segments(plotout[1,1,], y.axis, plotout[3,1,], y.axis, 
            lwd =  1.5, lty=2)#coef +/-1.96*se = 95% interval, lwd adjusts line thickness
points(plotout[2,1,], y.axis, pch = 19, cex=.6)
#emphsize winners
segments(winout[1,1,], y.axis, winout[3,1,], y.axis, lwd =  1.5)#coef +/-1.96*se = 95% interval, lwd adjusts line thickness
points(winout[2,1,], y.axis, pch = 19, cex=.7)
#axis for drawing axis
axis(1, at = c(low,hgh), tick=T, labels=c("",""),lwd.ticks=0, cex.axis = .8, mgp = c(2,.5,0))
axis(1, at = c( -.5, -.25, 0,  .25, .5), tick = T, cex.axis = .8, mgp = c(2,.5,0))
axis(2, at = y.axis, label = modname, las = 1, tick = T, cex.axis = .9) 
lines (c(0,0), c(0,max(y.axis)+.2), lwd=1, lty=3)
legend(-.65, 12,c("Lead donor weights"), bty="n", text.font=2)
legend(-.65, 8,c("Even Weights"), bty="n", text.font=2)
legend(-.65, 4, c("Geographic weights"), bty="n", text.font=2)
legend(-.4, 11, c("Best fitting models", "Loser in Clarke test"), lty=c(1,2), lwd=1,
            cex=.8)
dev.off()


#Substantive effects

#do 
#all channels without lead donor
#all channels with lead donor

sarhat <- function(pars, x, w){
    rho <- pars[1]
    k <- ncol(x)
    b <- as.matrix(pars[3:(k+2)])
    nd <- nrow(w)
    n <- nrow(x)
    d <- n/nd
    yhat <- rep(NA, n)
    
    for(i in 1:d){
        Ainv <- solve(diag(1, nd, nd)-rho*w[,,i])
        yhat[(1+nd*(i-1)):(nd+nd*(i-1))] <- Ainv %*% as.matrix(x[((i-1)*nd+1):((i-1)*nd+nd),])%*%b
    }
    return(yhat)
}

nosarhat <- function(pars, x){
    k <- ncol(x)
    b <- as.matrix(pars[3:(k+2)])
    yhat<- as.matrix(x)%*%b
}

###########################
#Total aid without lead donor
setwd("C:/Users/Martin Steinwand/Documents/Project Lead Donor/Data")

rawdat <- read.dta("DeliveryDyads.4.2.dta")
#ensure data is sorted on r_ccode year d_ccode:
if (rawdat$r_ccode[1]!=rawdat$r_ccode[2]) cat("\n\n WARNING, data not sorted properly!!!\n\n")
if (rawdat$d_ccode[1]==rawdat$d_ccode[2]) cat("\n\n WARNING, data not sorted properly!!!\n\n")
if (rawdat$year[1]!=rawdat$year[2]) cat("\n\n WARNING, data not sorted properly!!!\n\n")

########################
#Weak lead donors & no lead donors
workdat <- data.frame(rawdat[!rawdat$d_ccode==732 #& (rawdat$clead5of9r2==1| rawdat$clead5of9r2==0)
                & rawdat$ldstrong==0
                & is.na(rawdat$ODAdis_f)==0 & is.na(rawdat$population)==0
                & is.na(rawdat$voteshare_f)==0 & is.na(rawdat$rgdpch)==0,] )
attach(workdat)

vars <- cbind(ODAdis_f, 1, ODADonYear*10e-2, hhi, OilExp_rf*1e-8, 
TotImp_rf*1e-8, Camerica,Samerica, Africa, Meast, Asia, 
Oceania, population*1e-8, rgdpch*1e-2, colown, voteshare_f)

varspop <- cbind(AllPop, 1, ODADonYear*10e-2, hhi, OilExp_rf*1e-8, 
TotImp_rf*1e-8, Camerica,Samerica, Africa, Meast, Asia, 
Oceania, population*1e-8, rgdpch*1e-2, colown, voteshare_f)

#remove missing country dummies
vnames <- c("Aid", "Constant", "Total Donor Aid ", "Donor Concentration", 
                        "Oil Exports", "Total Imports", "Central America",
                         "South America", "Africa", 
                        "Middle East", "Asia","Oceania","Population", 
                        "GDP per capita, ppp","Former Colony", 
                        "Joint UN Voting")
colnames(vars) <- vnames

#get aggregates for 
nvars <- aggregate(vars, list(d_ccode, r_ccode), mean)
nvarspop <- aggregate(varspop, list(d_ccode, r_ccode), mean)
colnames(nvars)[1:2] <- c("d_ccode","r_ccode")
colnames(nvars)[ncol(nvars)] <- c("Joint UN Voting")
y <- nvars[,3]
ypop <- nvarspop[,3]
x <- nvars[,4:ncol(nvars)]

#counts 
n <- nrow(nvars)
nd <- length(table(nvars$d_ccode))
d <- n/nd


#W matrix -- even
#array
wall.ev <- array(NA, c(nd, nd, d))
wraw <- matrix(c(rep(1/(nd-1),nd)),nd, nd, byrow=TRUE)
diag(wraw) <- 0
wall.ev[,,1:d] <- wraw

yhat1 <- sarhat(res[[1]]$par, x, wall.ev)
yhatno1 <- nosarhat(res[[1]]$par, x)
dyhat1 <- yhat1 -yhatno1

###############################
#Total aid with strong lead donor
detach(workdat)
rm(workdat)

workdat <- data.frame(rawdat[!rawdat$d_ccode==732 & rawdat$ldstrong==1  & is.na(rawdat$ldstrong)==0
#                        & rawdat$nogovaid==0
                        & is.na(rawdat$ODAdis_f)==0 & is.na(rawdat$population)==0
                        & is.na(rawdat$voteshare_f)==0 & is.na(rawdat$rgdpch)==0,])
attach(workdat)

vars <- cbind(ODAdis_f, 1,  stronglead5of9r2, ODADonYear*10e-2, hhi, OilExp_rf*1e-8, 
        TotImp_rf*1e-8, Camerica,Samerica, Africa, Meast, Asia, 
        Oceania, population*1e-8, rgdpch*1e-2, colown, voteshare_f)

varspop <- cbind(AllPop, 1,  stronglead5of9r2, ODADonYear*10e-2, hhi, OilExp_rf*1e-8, 
        TotImp_rf*1e-8, Camerica,Samerica, Africa, Meast, Asia, 
        Oceania, population*1e-8, rgdpch*1e-2, colown, voteshare_f)

#remove missing country dummies
vnames <- c("Aid", "Constant", "Lead Donor", "Total Donor Aid ", "Donor Concentration", 
                                "Oil Exports", "Total Imports", "Central America",
                                 "South America", "Africa", 
                                "Middle East", "Asia","Oceania","Population", 
                                "GDP per capita, ppp","Former Colony", 
                                "Joint UN Voting")
colnames(vars) <- vnames

#get aggregates for 
nvars <- aggregate(vars, list(d_ccode, r_ccode, leadid5of9r2), mean)
nvarspop <- aggregate(varspop, list(d_ccode, r_ccode, leadid5of9r2), mean)
colnames(nvars)[1:3] <- c("d_ccode","r_ccode","leadid5of9r2")
colnames(nvars)[ncol(nvars)] <- c("Joint UN Voting")
y <- nvars[,4]
ypop <- nvarspop[,4]
x <- nvars[,5:ncol(nvars)]

#counts 
n <- nrow(nvars)
nd <- length(table(nvars$d_ccode))
d <- n/nd


#W matrix -- lead donor

#blog diagonal matrix or array
#get aggregated donor list 
dlist <- nvars[1:nd,1]

wall.ld <- array(NA, c(nd, nd, d))
wraw <- matrix(c(rep(.5/(nd-2),nd)),nd, nd, byrow=TRUE)
ldrow <- rep(1/(nd-1), nd)
for (i in 1:d){
    wsub <- wraw
    wsub[dlist==nvars[(1+(i-1)*nd),3],] <- ldrow
    wsub[,dlist==nvars[(1+(i-1)*nd),3]] <- rep(.5, nd)
    diag(wsub) <- 0
    wall.ld[,,i] <- wsub
}    
 

yhat2 <- sarhat(res[[3]]$par, x, wall.ev)
yhatno2 <- nosarhat(res[[3]]$par, x)
dyhat2 <- yhat2 -yhatno2


#Now plot effects 


bars <- c(median(dyhat1),median(dyhat2))*1000 #1000$ per 1000 people

postscript(file="C:/Users/Martin Steinwand/Documents/Project Lead Donor/Analysis/R&R.IO/SubstanceEffect.R&R.eps",
        onefile=FALSE, horizontal=FALSE, paper="letter", height=6, width=6)
barplot(bars,     space=.6, ylab="Change in allocation, US Dollars per 1000 capita", 
        ylim=c(-200,300), names.arg=c("All channels\n No Lead Donor", 
            "All channels\n Lead Donor"), cex.names=1.2 )
dev.off()


postscript(file="C:/Users/Martin Steinwand/Documents/Project Lead Donor/Analysis/FinalPaper/EffectSize.eps",
        onefile=FALSE, horizontal=FALSE, paper="letter", height=4, width=8)
barplot(t(bars), width=rep(.8,8),density=c(5, 5, 5, 5, 20, 20, 20, 20), 
        angle=c(45, 0, -45, 90, 45, 0, -45, 90), ylab="Share of Strong Lead Donors", 
        ylim=c(0,1), cex.names=.8, cex.lab=1.2, xlim=c(0,8), names.arg=c("Central \nAmerica", 
        "South \nAmerica", "Europe", "Africa", "Middle \nEast", "Asia", "Oceania"))
legend.lb(6.7,1, legend=dname, density=c(5, 5, 5, 5, 20, 20, 20, 20), 
        angle=c(45, 0, -45, 90, 45, 0, -45, 90), fill=TRUE, bty="n")
dev.off()



#graphical display of other results 
colnames(bds) <- c("rho", "var", "Constant", "Total Donor Aid ", "Donor Concentration", 
                                "Oil Exports", "Total Imports",  
                                "Central America", "South America", "Africa", 
                                "Middle East", "Asia","Oceania", "Population", 
                                "GDP per capita, ppp","Former Colony", 
                                "Joint UN Voting")



####################################
### Control variables

##Control vars total aid
##
modname <- c("No lead donor",
                "Lead donor")

modord <- c(1, 3)

plotout <- array(NA, c(3, ncol(bds[[modord[2]]])-3, length(modord)))
for (j in 1:length(modord)){
    if(ncol(bds[[modord[j]]])<ncol(bds[[modord[2]]])){
        plotout[,2:(ncol(bds[[modord[2]]])-3),j] <- unlist(bds[modord[j]])[10:length(unlist(bds[modord[j]]))]
    } else plotout[,,j] <- unlist(bds[modord[j]])[10:length(unlist(bds[modord[j]]))]
}


y.axis <- c(2:1)
vartit <- c("Lead donor dummy", "Total aid per donor & year", "Donor concentration",
            "Oil exports to donor", "Imports from donor", "Population", "GDP p.c.",
            "Former colony", "UN Voting")

low <- min(plotout[,4:7,], plotout[,14:17,], na.rm = TRUE)
hgh <- max(plotout[,4:7,], plotout[,14:17,], na.rm = TRUE)

lowall <- rep(NA, length(vartit))
hghall <- rep(NA, length(vartit))

varn <- c(4:7, 14:17)
varn <- c(1:5, 12:15)
for (k in 1:length(vartit)){
    lowall[k] <- min(plotout[,varn[k],], na.rm=TRUE)
    hghall[k] <- max(plotout[,varn[k],], na.rm=TRUE)
}

#save
setwd("C:/Users/Martin Steinwand/Documents/Project Lead Donor/Analysis/R&R.IO")
postscript(file="VarsAllAid.R&R.eps",onefile=FALSE,
 horizontal=FALSE, paper="letter", height=6, width=9)
par(mfrow=c(ceiling(sqrt(length(varn))), ceiling(sqrt(length(varn)))))
for(j in 1:length(vartit)){
    if(j==1|j==4 | j==7) par(mar=c(4, 8,2, .5), xpd=TRUE)else par(mar=c(4, 4, 2, 2.75), xpd=TRUE)
    plot(c(lowall[j], hghall[j]), c(1,max(y.axis)), type = "n", 
        axes=F,
        xlab = "", 
        ylab = "", #turn off axes and labels. 
            xlim=c(lowall[j], hghall[j]), c(1, max(y.axis)), xaxs = "r",
             main = vartit[j] ) 
    #the 3 lines below create horiztonal lines for 95% confidence intervals, and vertical ticks for 90% intervals
    segments(plotout[1,varn[j],], y.axis, plotout[3,varn[j],], y.axis, 
                lwd =  1.5)#coef +/-1.96*se = 95% interval, lwd adjusts line thickness
    points(plotout[2,varn[j],], y.axis, pch = 19, cex=.6)
    axis(1, at = c(lowall[j], hghall[j]), tick=T, labels=c("",""),lwd.ticks=0, 
        cex.axis = .8, mgp = c(2,.5,0))
    axis(1, #at = round(seq(lowall[j], hghall[j], length.out=5),1), 
        tick = T,cex.axis = .8, mgp = c(2,.5,0))
    if (j==1|j==4 | j==7)axis(2, at = y.axis, label = modname, las = 1, tick = T, cex.axis = .9) 
    if (lowall[j]<0) lines (c(0,0), c(.5,max(y.axis)+.2), lwd=1, lty=3)
    if (j==1|j==4 | j==7){par(font=2)
        legend(-130, 12,c("Total aid"), bty="n")
        legend(-130,8,c("NGO Aid"), bty="n")
        legend(-130, 4, c("Government Aid"), bty="n")
        par(font=1)
    }
}
dev.off()
